Non-Markovian Accumulators, Transition States, and Backwards Conversion
params =
list(
t_names = c("natural_history"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Healthy", "CVD", "Dead"), # State names
n_states = 3, # Number of states
n_cohort = 1, # Cohort size
n_cycles = 100, # Number of cycles in model.
cycle = 1, # Cycle length
initial_age = 55, # Cohort starting age
r_H_CVD = 0.15, # Rate of healthy -> CVD
hr_CVD = 10, # Hazard Ratio: CVD Death
r_H_D = 0.01 # Rate of healthy -> dead
)params$mR <-
with(params,{
r_CVD_D <- hr_CVD * r_H_D
R_ <-
array(data = c(0, 0, 0,
r_H_CVD,0, 0,
r_H_D,r_H_D + r_CVD_D,0,
r_H_D, 0,0),
dim = c(n_states, n_states, n_treatments),
dimnames = list(from = s_names,
to = s_names,
t_names))
R <- apply(R_,3, simplify=FALSE,function(x) {
diag(x) = -rowSums(x)
x * cycle
})
R
})Augmenting the transition matrix with non-markovian elements can address each question:
| Healthy | CVD | Dead | |
|---|---|---|---|
| Healthy | 0.852 | 0.131 | 0.017 |
| CVD | 0.000 | 0.896 | 0.104 |
| Dead | 0.000 | 0.000 | 1.000 |
| Healthy | CVD | Dead | trCVD | accCVD | |
|---|---|---|---|---|---|
| Healthy | 0.852 | 0.131 | 0.017 | . | . |
| CVD | 0 | 0.896 | 0.104 | . | . |
| Dead | 0 | 0 | 1 | . | . |
| trCVD | . | . | . | . | . |
| accCVD | . | . | . | . | . |
| Healthy | CVD | Dead | |
|---|---|---|---|
| Healthy | -0.16 | 0.15 | 0.01 |
| CVD | 0.00 | -0.11 | 0.11 |
| Dead | 0.00 | 0.00 | 0.00 |
We’ll start with the original transition rate matrix:
Next, we will add both a column and a row for the accumulator that tracks movement into the CVD health state.
For now, just include zeros.
| Healthy | CVD | Dead | accCVD | |
|---|---|---|---|---|
| Healthy | 0.852 | 0.131 | 0.017 | 0.139 |
| CVD | 0 | 0.896 | 0.104 | 0 |
| Dead | 0 | 0 | 1 | 0 |
| accCVD | 0 | 0 | 0 | 1 |
| cycle | Healthy | CVD | Dead | accCVD |
|---|---|---|---|---|
| 0 | 100000 | 0 | 0.0 | 0 |
| 1 | 85214 | 13107 | 1678.5 | 13862 |
| cycle | Healthy | CVD | Dead | accCVD |
|---|---|---|---|---|
| 0 | 100000 | 0 | 0.0 | 0 |
| 1 | 85214 | 13107 | 1678.5 | 13862 |
CVD CVD Death as well as death from background causes.accCVD accCVD transition probability.We’ll start with the original transition rate matrix:
| Healthy | CVD | Dead | trCVDDeath | |
|---|---|---|---|---|
| Healthy | 0.852 | 0.131 | 0.017 | 0.007 |
| CVD | 0.000 | 0.896 | 0.104 | 0.095 |
| Dead | 0.000 | 0.000 | 1.000 | 0.000 |
| trCVDDeath | 0.000 | 0.000 | 0.000 | 0.000 |
Healthy state (compound transitions) and the CVD state!| cycle | Healthy | CVD | Dead | trCVDDeath |
|---|---|---|---|---|
| 0 | 100000 | 0 | 0.0 | 0.00 |
| 1 | 85214 | 13107 | 1678.5 | 685.83 |
trCVDDeath and apply as usual to calculate accumulated costs from a one-time CVD death cost.# Define the cost vector
cost_payoff = c("Healthy" = 0, "CVD" = 500, "Dead" = 0, "trCVDDeath" = 2000)
# Sum up total costs in each cycle
costs_cycle = trace[["natural_history"]][-1,] %*% cost_payoff
# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8), 49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_costs = cycle_adj(costs_cycle)
total_costs[1] 5926.6
We’ll start with the original transition rate matrix:
$natural_history
to
from Healthy CVD Dead
Healthy -0.16 0.15 0.01
CVD 0.00 -0.11 0.11
Dead 0.00 0.00 0.00
R <-
cbind(R_,"tunCVDy1" = c(0,0,0), "tunCVDy2" = c(0,0,0),"NULL" = c(0,0,0)) %>%
rbind(.,"tunCVDy1" = c(0,0,0,0,0,0), "tunCVDy2" = c(0,0,0,0,0,0), "NULL" = c(0,0,0,0,0,0))
R Healthy CVD Dead tunCVDy1 tunCVDy2 NULL
Healthy -0.16 0.15 0.01 0 0 0
CVD 0.00 -0.11 0.11 0 0 0
Dead 0.00 0.00 0.00 0 0 0
tunCVDy1 0.00 0.00 0.00 0 0 0
tunCVDy2 0.00 0.00 0.00 0 0 0
NULL 0.00 0.00 0.00 0 0 0
P_ = expm(R)
P_["tunCVDy1","NULL"] = P_["CVD","Dead"]
P_["tunCVDy1","tunCVDy1"] = 0
P_["tunCVDy1","tunCVDy2"] = 1 - P_["tunCVDy1","NULL"]
P_["tunCVDy2","tunCVDy2"] = 0
P_["tunCVDy2","NULL"] = 1
P = P_[-which(rownames(P_)=="NULL"),-which(colnames(P_)=="NULL")]
P Healthy CVD Dead tunCVDy1 tunCVDy2
Healthy 0.85214 0.13107 0.016785 0.13862 0.00000
CVD 0.00000 0.89583 0.104166 0.00000 0.00000
Dead 0.00000 0.00000 1.000000 0.00000 0.00000
tunCVDy1 0.00000 0.00000 0.000000 0.00000 0.89583
tunCVDy2 0.00000 0.00000 0.000000 0.00000 0.00000
| Healthy | CVD | Dead | tunCVDy1 | tunCVDy2 | |
|---|---|---|---|---|---|
| Healthy | 0.852 | 0.131 | 0.017 | 0.139 | 0.000 |
| CVD | 0.000 | 0.896 | 0.104 | 0.000 | 0.000 |
| Dead | 0.000 | 0.000 | 1.000 | 0.000 | 0.000 |
| tunCVDy1 | 0.000 | 0.000 | 0.000 | 0.000 | 0.896 |
| tunCVDy2 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
A B C D
A 0.721 0.202 0.067 0.010
B 0.000 0.581 0.407 0.012
C 0.000 0.000 0.750 0.250
D 0.000 0.000 0.000 1.000
V <- eigen(m_P)$vectors
iV <- solve(V)
Ap <- iV %*% m_P %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
rownames(R) <- c("A", "B", "C", "D")
colnames(R) <- c("A", "B", "C", "D")
round(R,3) A B C D
A -0.327 0.311 0.002 0.013
B 0.000 -0.543 0.615 -0.072
C 0.000 0.000 -0.288 0.288
D 0.000 0.000 0.000 0.000
A B C D
A -0.327 0.311 0.002 0.013
B 0.000 -0.543 0.615 -0.072
C 0.000 0.000 -0.288 0.288
D 0.000 0.000 0.000 0.000
expm::logm Higham method returns same as Eigenvalue method. to
from 1 2 3 4
1 721 202 67 10
Call:
msm(formula = state ~ years, subject = ptnum, data = x, qmatrix = Q.init)
Maximum likelihood estimates
Transition intensities
Baseline
State 1 - State 1 -0.3271 (-0.3680,-0.2907)
State 1 - State 2 0.3271 ( 0.2907, 0.3680)
State 2 - State 2 -0.6454 (-0.8181,-0.5092)
State 2 - State 3 0.6454 ( 0.5092, 0.8181)
State 3 - State 3 -0.3982 (-0.7562,-0.2097)
State 3 - State 4 0.3982 ( 0.2097, 0.7562)
-2 * log-likelihood: 1572.2
A B C D
A 0.721 0.202 0.067 0.010
B 0.000 0.581 0.407 0.012
C 0.000 0.000 0.750 0.250
D 0.000 0.000 0.000 1.000
State 1 State 2 State 3 State 4
State 1 0.721 0.202 0.067 0.010
State 2 0.000 0.524 0.384 0.092
State 3 0.000 0.000 0.672 0.328
State 4 0.000 0.000 0.000 1.000
A B C D
[1,] 721 202 67 10
State 1 State 2 State 3 State 4
[1,] 721 202 67 10
log numerically unstable converting from probability to rate.log methods (eigen/Higham) usually results in impossible models.msm.